clear
clc

counter_maxs_cost=zeros(68,6);
counter_subs_cost=zeros(68,6);
counter_mergers_cost=zeros(68,6);
counter_maxs_comp=zeros(68,5);
counter_subs_comp=zeros(68,5);
counter_mergers_comp=zeros(68,5);

parfor zz=1:68
tic    
%%%%%loading data    
w=load('competitions.mat');
COMP=w.comp_list(zz,1);
x=load('data.mat');
y=load('qestimates.mat');
competitionid=x.competitionid;
maxscore = x.maxscore(competitionid==COMP,1);
pubscore_normal = x.pubscore_normal(competitionid==COMP,1);
length_days= x.length_days(competitionid==COMP,1);
length_days=length_days(1);
rewardquantity= x.rewardquantity(competitionid==COMP,1);
rewardquantity=rewardquantity(1);
mergers_data= x.totmergers(competitionid==COMP,1);
mergers_data=mergers_data(1);
nSim = 500;

%%%%%state space parameters
N=40; %number of individuals
T=360; %number of periods
Svalues=unique([unique(maxscore);max(maxscore)+[0.002:0.002:0.08]']);
nScore = 20;
Svalues = Svalues(3:nScore+2,1);
nsubs_data=size(maxscore,1);
lambda1=0.75;
lambda2=1-lambda1;

%%%%%q-function
betachangehat=[y.beta0(y.cid==COMP,1) ; y.beta1(y.cid==COMP,1); y.beta1tr(y.cid==COMP,1)];
probchange = @(score) (exp(betachangehat(1) + betachangehat(2)*score)./(1+exp(betachangehat(1) + betachangehat(2)*score)));
q_sp=probchange(Svalues);
%q_sp(end,1)=0;
probchange = @(score) (exp(betachangehat(1) + betachangehat(3) + betachangehat(2)*score)./(1+exp(betachangehat(1)+ betachangehat(3) + betachangehat(2)*score)));
q_team=probchange(Svalues);
%q_team(end,1)=0;

%load estimates
asdf=load(sprintf('Estimates/%s_%02d.mat', 'estimates',COMP));
sigmaS=asdf.sigmaS;
sigmaT=asdf.sigmaT;

%Counterfactual cost of team formation
sigma_tilde = @(sigma,alpha) (alpha*sigma/(1+sigma))/(1-(alpha*sigma/(1+sigma)));
for ww=1:6
alph = 0.5+(ww-1)*0.1;    
[mergers,subs,maxs,~]=equilibrium(sigmaS,sigma_tilde(sigmaT,alph),N,T,lambda1,lambda2,q_team,q_sp,nScore,Svalues,mergers_data,nSim);
counter_maxs_cost(zz,ww) = maxs;
counter_mergers_cost(zz,ww) = mergers;
counter_subs_cost(zz,ww) = subs;
end

%Counterfactual competition
for ww=1:1:5
    NN=36+2*(ww-1);
    TT=324+18*(ww-1);
[mergers,subs,maxs,~]=equilibrium(sigmaS,sigmaT,NN,TT,lambda1,lambda2,q_team,q_sp,nScore,Svalues,mergers_data,nSim);
counter_maxs_comp(zz,ww) = maxs;
counter_mergers_comp(zz,ww) = mergers;
counter_subs_comp(zz,ww) = subs;
end
end

writematrix([counter_mergers_cost,counter_subs_cost,counter_maxs_cost],'counter_out_cost.csv') 
writematrix([counter_mergers_comp,counter_subs_comp,counter_maxs_comp],'counter_out_N.csv') 

main_inference